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The collective dynamics of excitatory pulse coupled neurons with spike timing dependent plasticity 
(STDP) is studied. The introduction of STDP induces persistent irregular oscillations between 
strongly and weakly synchronized states, reminiscent of brain activity during slow-wave sleep. We 
explain the oscillations by a mechanism, the Sisyphus Effect, caused by a continuous feedback 
between the synaptic adjustments and the coherence in the neural firing. Due to this effect, the 
synaptic weights have oscillating equilibrium values, and this prevents the system from relaxing into 
a stationary macroscopic state. 
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Sisyphus was the mythological king of Corinth com- 
pelled to roll a heavy boulder up a hill, only to watch it 
roll back down as it approached the top. Sisyphus was 
condemned by Zeus for his iniquity and pride to repeat 
eternally his efforts, without any hope of success. How- 
ever, in the brain such endless motion can have a positive 
functional relevance. Fluctuating spontaneous activity 
has been observed in several areas of the brain (l|. In 
particular, irregular oscillations between more and less 
synchronized states have been revealed in the hippocam- 
pus during slow-wave sleep and this activity has been 
related to memory consolidation in the neocortex [2j. 

Recent studies have suggested synaptic plasticity as a 
fundamental ingredient to ensure multistability in neu- 
ronal circuits [3|, y, la, [6| ■ In particular, spike timing de- 
pendent plasticity (STDP) is considered one of the cen- 
tral mechanisms underlying information elaboration and 
learning in the brain [2|. A series of experiments per- 
formed in vivo and in vitro on neural tissues revealed 
that the strength of a synapse, conveying spikes from 
a presynaptic to a postsynaptic neuron, depends cru- 
cially on the precise spike timing of the two connected 
neurons [8l4l0|. The STDP rules prescribe that when- 
ever the presynaptic (postsynaptic) neuron fires before 
the postsynaptic (presynaptic) one, the synapse is po- 
tentiated (depressed). The synapse is modified only if 
the spikes occur within certain time intervals [learning 
windows). Asynrmetric learning windows have repeat- 
edly been found experimentally (e.g. see plTll3j ). This 
asymmetry is a prerequisite, at least in phase oscillator 
networks, to observe the coexistence of states character- 
ized by different levels of synchrony [3J, y] . Furthermore, 
in the presence of propagation delays STDP can provide 
a negative feedback mechanism contrasting highly syn- 
chronized network activity and promoting, in randomly 
driven networks, the emergence of states at the border 



between randomness and synchrony [5(. 

In this Letter a novel deterministic mechanism, the 
Sisyphus effect (SE), able to generate spontaneous fluc- 
tuations in a neural network between asynchronous and 
synchronous regimes is presented. In particular, we study 
excitatory pulse coupled neural networks with STDP, 
where the interaction among neurons is mediated by a- 
pulses [14| . For non plastic interactions, the excitatory 
coupling leads to synchronization only for sufficiently fast 
synapses [15j. Furthermore, the desynchronizing effect is 
amplified at large coupling |16j. In absence of plastic- 
ity the macroscopic activity of the network is station- 
ary: asynchronous for large synaptic weights and par- 
tially synchronized for sufficiently weak coupling [ly, |_17| . 

The introduction of STDP completely modifies the dy- 
namical landscape leading to a regime where a strongly 
and a weakly synchronized state coexist. The activity of 
the network is thus characterized by irregular oscillations 
between these two states. These transitions are driven by 
the evolution of the synaptic weights, which in turn is dic- 
tated by the level of synchrony in the network. For small 
synaptic weights the system is fully synchronized, while 
above a critical coupling it desynchronizes. Further- 
more, whenever the network is synchronized (desynchro- 
nized) the synaptic weights tend towards large (small) 
equilibrium values corresponding to asynchronous (syn- 
chronous) dynamics. In summary, the neuronal activity 
can be represented in terms of an order parameter dif- 
fusing over an effective free energy landscape displaying 
two coexisting equilibrium states. Small (large) synaptic 
weights tilt the landscape towards the strongly (weakly) 
synchronized state, in turn the induced activity increases 
(reduces) the weights until a tilt in the opposite direction 
occurs. Thus the landscape oscillates endlessly. 

The model. We study a fully coupled network of N 
Leaky Integrate-and-Fire neurons, for which the mem- 



brane potential Vi(t) €[0:1] of neuron i evolves as: 



Vi(t)=a-Vi(t) + Ii(t) 



1. 



,N 



(1) 



whenever the neuron reaches the threshold Vt = 1, an a- 
pulse p a {t) = oi 2 £exp(— at) is instantaneously transmit- 
ted to all other neurons and Vi is reset to zero. Further- 
more, a > 1 is the suprathreshold DC current, li = gEi 
the synaptic current, and g the excitatory homogeneous 
coupling. The field Ei represents the linear superposi- 
tion of the pulses received by neuron i and its evolution 
is ruled by a second order ODE (Eq. (S2) in |18j). For 
a fully coupled non-plastic network the synaptic weights 
associated to the connection from the presynaptic j'-th 
neuron to the postsynaptic i-th one are w%j = 1 (apart 
from the autaptic terms: wu = 0). 

In the presence of plasticity, we assume that the 
weights evolve in time accordingto a nearest-neighbour 
STDP rule with soft bounds [2|45j, Ig • Therefore in the 
case of a post- (presynaptic) spike, emitted by neuron i 
(j) at time t, the weight Wij is potentiated (depressed) 



as w, 



i(t + ) 



"', 



r«(*) = 4 



,-(£") +ry(£), with 

6 ij 

p[w M - Wij(t~)]e T + if <% > 
—dwij(t )e T - if Sij < 



(2) 



where 5y = t—t^ > (5 zj = t^-t < 0) is the firing time 
difference and t^ the last firing time of neuron k. The 
potentiation and depression factors (p and d, resp.) co- 
incide, unless otherwise specified [22j ]. The bounds keep 
the synapses from achieving unrealistically large values or 
becoming inhibitory, namely < Wji < Wm- The learn- 
ing windows over which post- (pre-) synaptic spikes will 
cause synaptic potentiation (depression) are indicated as 



11 



we as- 



r+ (t_ ). Following experimental evidences 
sume r_ > t+ . The degree of synchronization of the neu- 
rons is measured by the order parameter |9j, [l0| R{t) = 



N ^A 



i6k(t)\ 



where k {t) = 2 7 r(i-C J )/( i 



(k)s/,.(k) 



m+l 



t { * ] ) is 



the phase of the fc-th neuron at time t between its m-th 
and (m + l)-th spike emission. A perfectly synchronized 
(asynchronous) system has R = 1 (R = 0), while inter- 
mediate values indicate partial synchronization. 

Phase Diagram. We analyze how the phase diagram 
of the network is modified by the plasticity. In particu- 
lar, we focus on the variation of the neuronal coherence 
by varying the DC current. Similar results can be ob- 
tained by varying the coupling g and the pulse width a 
(as shown in [18||). To compare with previous results ob- 
tained without plasticity, we fix g = 0.4 and a = 9 as 
in [7|, [l7| . In absence of plasticity the homogeneous sys- 
tem exhibits two phases: an asynchronous regime with 
R = 0, and a partially synchronized phase with finite 
17j |. The emergence of one or the other regime de- 



R 



pends crucially on the ratio of two time scales: the pulse 
rise time 1/q and the interspike-interval (ISI) [7l.ll5|. For 



slow synapses (relative to the ISI) the system dynamics 
is asynchronous, while for sufficiently fast synapses co- 
herent oscillations emerge. The system becomes fully 
synchronized only for instantaneous synaptic rise times 
(i.e. a — > oo). For fixed network size N and pulse shape, 
the ISI can be reduced by increasing either the external 
DC current or the synaptic coupling. Therefore partial 
synchronization is observable for sufficiently small a or 
g values (whenever a > 3.4), while incrementing these 
parameters will desynchronize the system [161 ] (as shown 
in Fig. la in (li|). 




FIG. 1: (color online) (a) Average order parameter R as a 
function of a for the non plastic network (black solid line) 
and in the presence of STDP for TV = 200 (red filled squares) 
and TV = 500 (blue dashed line) (b) Free energy profile F(R) 
versus R for N = 200, obtained by evaluating R at regular 
time intervals At = 1 for a time span ~5x 10 6 . (c) Time 
evolution of R(t) (black solid line) and of W(t) (red dashed 
line) for N = 500. The dotted (blue) line is the W predicted 
via Eq. (fSJl. The data refer to a — 1.3, g — 0.4, a — 9, 
d = p = 0.01, r_ = 3r + = 0.3, and wm — 2, and are measured 
after a transient ~ 10 5 . 

The average level of synchronization R is reported in 
Fig. \Hn as a function of a for the non-plastic and plas- 
tic cases. In absence of plasticity the system is partially 
synchronized for low DC currents and asynchronous for 
a > a c ~ 1.35. The introduction of plasticity does not 
alter the scenario at small a- values, where the system is 
in a high synchronization (HS) regime. The main differ- 
ence is observable in the dynamics of R(t), which dis- 
plays irregular oscillations: the associated Fourier spec- 
trum resembles a Lorcntzian with a small subsidiary peak 
around period ~ 34 — 36. However, for sufficiently large 
currents, namely a > 1.5, the asynchronous regime is 
substituted by a state of low synchronization (LS) char- 
acterized by a rapidly fluctuating order parameter (over 
a time scale of the order of 70 - 150) with an associated 
small level of synchronization R ~ 0.32 ± 0.12. At inter- 
mediate a-values, in the range a £ [1-23; 1.46], R exhibits 
wide irregular temporal oscillations between values ~ 1 
and zero with characteristic time scales ~ 1100 — 1400. 
These latter oscillations represent low frequency fluctu- 
ations (LFFs), while rapid fluctuations are still present 



over time scales ~ 50 — 60 (see Fig. [TJ:). 

In this Letter we will mainly focus on the intermediate 
regime, fixing a = 1.30, where the LFF of R(t) resembles 
the evolution of a particle in a double well potential sub- 
ject to thermal fluctuations. To clarify this analogy we 
have estimated the PDF, P(R), of the order parameter, 
by examining its trajectory for a sufficiently long time 
span, and derived the associated free energy profile as 
F(R) = - logP(iJ). As shown in Fig. [TJj, F(R) exhibits 
two minima corresponding to a HS phase at Rh — 0.905 
and a LS state at Rl — 0.32. The 2 coexisting minima 
are separated by a saddle, located at Rs — 0.61. As 
clarified in the following the jumps between minima are 
driven by the macroscopic evolution of network plastic- 
ity. The rapid fluctuations, present in all regimes, are 
instead due to the microsopic evolution of the synaptic 
weights, which can be interpreted as a noise source for 
the dynamics of R(t) . The analysis of these noise-induced 
oscillations goes beyond the scope of this Letter and it is 
left for future studies. 

Constrained phase diagram As shown in Fig. [TJ; the 
LFFs of R(t) are associated with oscillations in the av- 
erage synaptic weight W(t) = J2i j Wji(t)/N(N — 1). In 
particular when the system is in a HS (LS) state W in- 
creases (decreases). 

To better investigate the origin of these correlations 
and the interaction between the STDP induced synap- 
tic dynamics and the level of synchronization in the sys- 
tem we perform the following numerical experiments. We 
simulate the system by constraining the synaptic weights 
to have a constant average value Wo, by rescaling, at reg- 
ular time intervals, the weights wtj. Initially, Wo = and 
we follow the evolution of the system for a time span T$. 
We then perform a new simulation for the same time 
lapse with a larger Wo value, starting from the last con- 
figuration of the previous run. The procedure is repeated 
by increasing Wo at regular steps AWo until W) = wm 
is reached. Then with the same protocol Wo is decreased 
(in steps of AWo) until finally Wo returns to zero [261 ] . 
The results of these simulations are shown in Fig. [2] for 
N = 200. At low W) the system is fully synchronized, 
while with increasing Wo the system desynchronizes via 
a discontinuous transition. By further increasing Wo the 
level of synchronization continues to decrease and an- 
other smooth transition seems to occur. For the expla- 
nation of the SE it is sufficient to limit the analysis to 
the first transition. 

As shown in Fig. [2 the constrained system exhibits a 
hysteretic transition from HS to LS (from LS to HS) for 



R 



W, 



(i) 



(2) 



0.76(5) (Wq = 0.65(5)) by increasing (decreas- 
ing) the control parameter Wq [27j. This implies that 



in the interval [W , W ] the two regimes coexist and 
that HS or LS is observable depending on the initial state 
of the network. 

Mean field synaptic evolution. In order to gain some 
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FIG. 2: (color online) R versus Wo as measured for increasing 
(black filled circles) and decreasing (empty red squares) Wo- 
The (blue) vertical dashed lines indicate the fixed point values 
Ws and Wa [laj. Results averaged over 8 different initial 
conditions, T s = 1,000, AWo = 0.02 (for clarity only one 
point every two is shown). (Inset) Conditional free energy 
profiles Fi(R) (black solid line) and Fd(R) (red dashed line) 
obtained during USs. Both curves are vertically shifted to 
achieve zero as minimal value. Parameters as in Fig. [l] and 
N = 200. 



insight into the evolution of the system during uncon- 
strained simulations (USs), let us consider a mean field 
equation for the synaptic weight evolution. The aver- 
age synaptic weight modification T, for each presynaptic 
spike, can be written as [5| 

T(t)=p{w M -W) dSP(S)e- r + -dW dSP(~S)c T - 



(3) 
where P(5) is the PDF of the time differences 6 between 
postsynaptic and presynaptic firing measured. To test 
the predictive value of Eq. ([3]), we have measured from 
an US P(S) at regular intervals At. By employing this 
information we can predict quite well the evolution of the 
synaptic weight as W(t+At) = W(t)+T(t) (sec Fig. Hh). 
By assuming that the postsynaptic neuron is firing 
with period Fq, we are able to derive the time difference 
distribution P(8) for the two limiting cases: fully syn- 
chronized and asynchronous dynamics. In the fully syn- 
chronized (asynchronous) situations we expect a distribu- 
tion of the form P s (S) =V(6)+V(6-T ) (P A (5) = 1/T ) 
defined in the interval [0 : To]. Here T> denotes a Dirac 
delta function. These guesses are essentially confirmed by 
direct USs as shown in Fig. 3 in [13]. Therefore in these 
two cases an analytical estimation of T can be obtained. 
Furthermore, in both cases Y vanishes for a finite value 
of the average synaptic weight, namely Ws (Wa) for the 
synchronized (asynchronous) situation. Furthermore, for 
W < Ws (W > Ws) the synapses are in average poten- 
tiated (depressed) and analogously for the asynchronous 
case. This implies that Ws (Wa) is a stable attractive 
point for the dynamics of W in the synchronized (asyn- 
chronous) regime (for a definition of Ws an d Wa see 



Eq. (S8) and (S10) in [lj). 

Sisyphus mechanism. We are now able to explain the 
behavior reported in Fig. QJ; for R(t) and W(t). Let us 
suppose that the system is in the HS phase with an as- 
sociated low coupling W < Wq . However, in this situa- 
tion the attractive fixed point Ws is above the transition 
point Wq (see Fig. [2]). Therefore W keeps increasing, 
until for W > W the system starts to desynchronize 
and to approach the LS state. In this phase the P(S) be- 
comes almost flat (see Fig. 3b in [18() and the attractive 
point for the synaptic evolution will be Wa, located be- 
low Wq . The motion towards Wa leads to a decrease of 
W . Whenever the average synaptic weight crosses W 
the neurons begins to resynchronize. Finally, the system 
will return to the HS state from where it started. The 
cycle will repeat indefinitely and is the essence of the SE. 

The above arguments arc approximate, because the 
system is never exactly fully synchronized or desynchro- 
nized, instead it passes through a continuum of states, 
each associated to a different fixed point in JU-spacc. 
The relevant aspect is that the fixed points associated 
to the HS (LS) phase are larger than the transition point 



(i) 



rVh 



(smaller than W ). As we have verified this is 



W i 

indeed the case, therefore the mechanism is still valid. 
To perform a direct test of the validity of our analysis, 
we have measured the PDF of R conditioned to the fact 
that W was increasing (decreasing) during an US. From 
these PDFs we derived the corresponding free energy pro- 
file Fj{R) (F D (R)), As shown in the inset of Fig. [2]^/ 
has a unique minimum at Rh , while Fjj has an absolute 
minimum at Rl and a shoulder around R ~ 0.8. These 
results confirm that the equilibrium attractive values for 
W are located opposite to the transition points, because 
when the system is in the HS (LS) regime the synaptic 
weights increase (decrease) continuously trying to reach 
the corresponding fixed points. 

The SE should be active whenever the transition val- 



-d) 



-12) 



ues Wq and Wq are both contained within the interval 
[Wa,Ws]- To verify this statement we have measured 



-d) 



(2) 



Wq , Wq and the fixed points for various DC currents 
within the interval < a < 2 (data shown in Fig. [3]). We 
observe that the transition is hysteretic in the interval 
a e]0; 1.40], while for larger values W 



i) 



' and Wq essen- 



tially coincide. Furthermore, W becomes larger than 
W s at a ~ 1.18, while W A > W {1) ,W {2) for a > 1.50. 
Thus we expect that F(R) exhibits two coexisting min- 
ima, due to the SE, when 1.18 < a < 1.50. To verify this 
conjecture we estimate the free energy barrier heights AF 
separating the HS and the LS state from the intermedi- 
ate saddle for various a values. As shown in the inset of 
Fig. [3l the barrier associated to the HS state diverges ex- 
ponentially when approaching a ~ 1.18. Therefore, the 
HS regime is the only possible at smaller a-valucs. On 
the other hand the two minima merge and the associated 
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FIG. 3: (color online) W versus the DC current a. The shaded 
area represents the VF-values measured during USs. The up- 
per black (lower red) dashed line refers to the estimated Wq 
(Wq ). The error bars have been evaluated over 5 differ- 
ent realizations of the constrained simulations. The upper 
(lower) solid blue line represents the fixed point values Ws 
(Wa)- (Inset) The red circles (black squares) refer to the 
free energy barrier AF separating the HS (LS) state from the 
saddle. The blue line is an exponential fit to the HS barrier 
height. The remaining parameters as in Fig. [l]and N = 500 



barriers vanish for a > 1.48 indicating that the LS state 
is the unique remaining at large a. Furthermore, the 
distributions of the W-values measured during USs are 
reported in Fig. [3] as a shaded area: these values include 
the transition interval [W (1) ; W (2) ] for 1.20 < a < 1.48. 

In conclusion, the SE should be observable in pulse 
coupled neural networks whenever the excitation has a 
desynchronizing effect. This is in general verified for any 
kind of neuronal response (type I or typ e II) for suffi- 
ciently slow synaptic interactions [la , Il6| . Furthermore, 
we have verified that the SE persists by setting p > d, as 



13] 



suggested by experimental evidences 
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ADDITIONAL MATERIAL 



THE MODEL 

We study a fully coupled network of N Leaky 
Intcgratc-and-Fire (LIF) neurons, for which the mem- 
brane potential Vi(t) €[0:1] of neuron i evolves as: 



where a > 1 is the suprathrcshold DC current, and /; 
the synaptic current due to the coupling with the rest 
of the network. Following [1[ , we assume that whenever 
the i-th neuron reaches the threshold Vi = 1, an a-pulse 
p a {t) = a 2 texp(—at) is instantaneously transmitted to 
all the other neurons and its membrane potential is re- 
set to Vi = 0. The synaptic current can be written as 
•A' CO = gEi(t), with g > representing the excitatory 
homogeneous coupling while the field Ei(t) is given by 
the linear superposition of the pulses received by neuron 
i in the past. For a-pulses the time evolution of the field 
Ei(t) is ruled by the following second order differential 
equation: 



Ei(t) + 2aEi{t) + a 2 Ei{t) 



N 



- Y^ Wij5(t-t„) 



n\t n <t 



. (5) 

with {t n } being the firing times until the present time 
t. For a fully coupled network, in absence of plasticity, 
the synaptic weights w%j appearing in Eq. ([5]) are all 
equal to one (apart from the autaptic terms which are 
set to zero). In presence of plasticity, we assume that 
the synaptic weights evolve in time according to a STDP 
rule with soft bounds, namely [2 4] 



Wij(t) =p{w r , 



w iJ {t))A j S i 



dwij(t)BiSj 



(6) 



where d (resp. p) is the potentiation (rcsp. depres- 
sion) amplitude, and S& represents the series of spik- 
ing times of neuron k until time t. The presence of the 
bounds implies that < Wij < w max . By assuming that 
the synapses have memory just of the last emitted spike 
[nearest-neighbour STDP rule), the variables Aj (resp. 
Bi) evolves as follows [2|, |4| 



r+Aj 



-Aj+O-A^Sj , r-Bi = -Bi+il-BJSi (7) 



where r + (resp. r_) are the time windows over which 
post- (resp. pre-) synaptic spikes will cause potentia- 
tion (resp. depression) of the synapse. As pointed out 
by Izhikevich and Desai [5| the nearest neighbour imple- 
mentation of the STDP rule is the only one consistent 
with the classical results on long-term potentiation and 
depression by Bicnenstock-Cooper-Munro [6[ 

Therefore in the case of a post-synaptic (pre-synaptic) 
spike, emitted by neuron i (j) at time t, the weight Wij 
is potentiated (depressed) as Wij(t + ) = Wij(t~) + Fij(t), 
with 

( _ _h 2 _ 

p[wm -Wij(t )]c T + if Sij > 

r«(t) = < (8) 



Vi(t)=a-Vi(t) + Ii(t) 



1. 



,N 



—dwij(t )c T - if 8^ < 

where 5^ = t - *Ci) > {5 rj = t® - t < 0) is the firing 
(4) time difference and t^ the last firing time of neuron k. 



Since the plasticity rule depends critically on the pre- 
cision of the spiking events, it is necessary to employ 
an accurate integration scheme to update the evolution 
equations. For this purpose we adopted an event-driven 
algorithm with event queue conjugating high accuracy 
with a fast implementation [7J, |8[ . 

The degree of synchronization of the neuronal popula- 
tion can be characterized in terms of the order parame- 
ter n [13 



R(t) = 



F^ 1 



where 



i$ k {t) 



(fch 



(9) 



h(t) = 2tt 



(t~t\rV) 



(t 



(fe) 
m+1 



t (k) ) 



(10) 



is the phase of the fc-th neuron at time t between its m- 
th and m + 1-th spike emission, occurring at times tin 
and t m+1 , respectively. A non zero R value is an indica- 
tion of partial synchronization, perfect synchronization is 
achieved for R = l, while a vanishingly small R ~ 1/yN 
is observable for asynchronous states in finite systems. 
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We report in Fig. |4]thc phase diagram for the non plas- 
tic fully connected network for N = 100 in the planes 
(g,a) and (g,a). In this case the asynchronous regime 
corresponds to a so-called splay state [l|, |j| : this is an 
exact solution for the system which is perfectly asyn- 
chronous (R = 0) . For this periodic solution we are able 
to perform an analytical stability analysis [7[ , reported 
in Fig. |U The (green) circles in Fig. |4] indicates the 
marginal stability line for the the splay state, it is known 
that whenever the splays state looses its stability it gives 
rise to a partially synchronized regime via a Hopf super- 
critical bifurcation [l[ . With reference to the parameter 
values considered in this Letter, partial synchronization 
emerges for a < a c ~ 1.35 for fixed coupling g = 0.4, and 
for g < g c ~ 0.4676 for fixed pulse width a = 9. 

In order to characterize the dynamical phases observed 
in the plastic and non plastic networks, we have also es- 
timated the order parameter R averaged over a certain 
time span. 

The behaviour of R as a function of g and a are re- 
ported in in Fig. [SJ similarly to the non plastic case the 
system desynchronizes also in presence of STDP for in- 
creasing synaptic coupling g and pulse rise times l/a. 
However, the perfectly asynchronous regime is substi- 
tuted by a low synchronization phase where R ~ 0.2. 
Furthermore in proximity of the transition region from 
high synchronous to the low synchronous regime R(t) 



FIG. 4: (color online) Phase diagram for the homogeneous 
model without plasticity in (a) we report the phase plane 
(g, a), while in (b) the phases are shown in the (g, a)-plane. 
The (green) filled circles indicate the critical values for which 
the asynchronous state (namely, the splay state) becomes un- 
stable; the dashed (magenta) line refers to g = 0.4 in (a) and 
a = 9 in (b). The other parameters are fixed to a = 9 in (a) 
and g — 0.4 in (b) and the system size is N = 100. 



exhibits large low frequency fluctuations as in the case 
studied in the Letter. 



TIME DIFFERENCE DISTRIBUTIONS 

The probability density distributions (PDF) for the time 
differences 5 between postsynaptic and presynaptic fir- 
ing time can be easily derived in two limiting case by 
assuming a constant inter-spike interval (ISI) To for the 
post-synaptic neuron. For fully synchronized neurons, 
we expect the pre- and postsynaptic neurons to fire to- 
gether S = or with a delay given by the ISI, namely 
5 = T . Therefore the distribution will be Ps(S) = 
V{5) + V{8 — T ), where V(x — xo) denotes a Dirac delta 
function centered in xq. The other situation we consider 
is that corresponding to perfect asynchrony, in this case 





FIG. 6: (color online) Probability distribution functions P(S) 
as obtained during an unconstrained simulation by consid- 
ering the TV x (N — 1) Sij values associated to the last N 
spikes preceeding a strongly (resp. weakly) synchronized state 
corresponding to an order parameter value R ~ 0.98 (resp. 
R ~ 0.11). The estimated PDF are reported in panel (a) and 
(b), respectively. The red dashed lines refer to Ps in (a) and 
Pa in (b) defined in the text. The parameters of the simula- 
tion are a = 9, r_ = 3r+ = 0.30, d — p = 0.01, g = 0.4 and 
the network size is JV = 500 



Asynchronous dynamics 



FIG. 5: (color online) (a) Average order parameter R versus 
a for the systems without plasticity (black circles) and with 
STDP (red squares); (b) R versus g for the systems without 
plasticity (black circles) and with STDP (blue squares). The 
parameters are a — 1.3 and N = 200, in (a) g — 0.4 and in 
(b) a = 9. The average have been performed over time spans 
~ 10 4 — 10 5 time units, after discarding transients of duration 
~ 10 4 -10 5 . 



we expect that 8 will take all values in the interval [0 : To] 
and all the values will be equiprobable therefore we ex- 
pcct P A {8) = l/T . 

In Fig. [6] we compare the predictions with the mea- 
sured distributions, the agreement is reasonable in view 
of the fact that the two considered states do not corre- 
spond exactly to R — 1 and R = and that the ISI is 
not constant over all the neuronal population. Therefore, 
we can consider Pa(8) and Ps(8) as reasonable approxi- 
mation of the true distributions in the two extreme cases 
achieved by the system during its evolution. 



In this situation P(S) = Pa{8) and we can rewrite Eq. 
(3) as follows 



T A 



To 



(w M -W)(l-e T +)-3VF(l-c --) 



(11) 

where we have assumed p = d and r_ = 3t + as for the 
most part of the simulations studied in the Letter. The 
quantity Ta vanishes for W = Wa and it is positive (resp. 
negative) for W < Wa (resp. W > Wa), therefore for 
the dynamics of W(t) 



W A = 



.In 



w M 1 - e T + 



To 

4-e T + 



(12) 



3e 



is a stable fixed point. 



Fully synchronized dynamics 



MEAN FIELD SYNAPTIC EVOLUTION 

By assuming that the post-synaptic neuron fires with 
constant period To we can perform the integrals appear- 
ing in Eq. (3) in the letter in the 2 limiting cases dis- 
cussed above. 



For the fully synchronized situation P(S) = Ps(S) and 
Eq. (3) becomes 



Ts =p 



w M (l+c T +)-W{2 + c T + +e T - ) 



(13) 



where we have assumed once more p = d and t_ = 3t + . 
The quantity T$ vanishes for W = Ws and it is positive 



8 



(resp. negative) for W < Ws (rcsp. W > Ws), therefore 
the solution 



W S 



w M I 1 + e T + 



_Zq_ _Zq_ 
2 + c T + +e T - 



represents an attractive fixed point. 
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